import numpy as np
import math

R = 8.314

def antoine(A, B, C, T):
    """计算安托因方程"""
    P = A - B / ( T + C )
    P = 10 ** P
    return P

def G_ex_MATRIX(x, tau, alpha, T):
    """计算过剩吉布斯自由能"""
    G = np.exp( -alpha * tau )
    G_ex = tau * G @ x  / ( G @ x ) * x
    G_ex = np.sum( G_ex, axis = 0 ) * R * T
    return G_ex

def gamma_cal(x1, y1, P, Pvs):
    """计算活度系数"""
    gamma = np.divide(P * y1 / Pvs, x1, out = np.zeros_like(x1), where=x1!=0)
    return gamma

def NRTL(x1, x2, tau12, tau21, alpha):
    """NRTL方程计算活度"""
    G12 = np.exp(-alpha * tau12)
    G21 = np.exp(-alpha * tau21)
    gamma = np.exp (x2 ** 2 * (tau21 * G21 ** 2 / (x1 + x2 * G21) ** 2 +\
                               tau12 * G12 / (x2 + x1 * G12) ** 2))
    return gamma

def NRTL_MATRIX(x, tau, alpha):
    """NRTL方程计算活度"""
    G = np.exp(-alpha * tau)
    term0 = tau * G @ x / ( G @ x ) #计算第一项，维度n * 1
    term1 = G * x #维度n*n，利用广播性质
    term2 = np.tile(G @ x, (1, np.size(x))) #在列上复制n次，维度n*n
    term3 = np.tile(term0, (1, np.size(x))) #在列上复制n次，维度n*n
    gamma = np.exp(term0.T + np.sum(term1 / term2 * ( tau - term3 ), axis = 0)) #按列求和
    return gamma


def HBT(Tr, Vc, w):
    "计算液体饱和蒸汽压状态下的密度"
    a = -1.52816
    b = 1.43907
    c = -0.81446
    d = 0.190454
    e = -0.296123
    f = 0.386914
    g = -0.0427258
    h = -0.0480645
    if Tr > 0.25 and Tr < 0.95:
        V0 = 1 + a * (1 - Tr) ** 0.333 + b * (1 - Tr) ** 0.667 +\
        c * (1 - Tr) + d * (1 - Tr) ** 1.333
    else:
        print(f"警告，计算VSL(V0)时Tr超出范围")
    if Tr > 0.25 and Tr < 1:
        Vd = (e + f * Tr + g * Tr ** 2 + h * Tr ** 3) / (Tr - 1.00001)
    else:
        print(f"警告，计算VSL(Vd)时Tr超出范围")
    VS = Vc * V0 * (1 - w * Vd)
    return VS
    
class CUBIC_EOS:
    """定义立方形状态方程中的各个变量"""

    def __init__(self, Pc, Tc, Z, w):
        "初始化常用参数"
        self.Pc = Pc
        self.Tc = Tc
        self.Z = Z
        self.w = w

    def PR_a(self, T):
        "计算PR方程的常数a"
        Tr = T / self.Tc
        F = 0.37464 + 1.54226 * self.w - 0.26992 * self.w ** 2
        a = (1 + F * (1 - Tr ** 0.5)) ** 2
        a = 0.457235 * R ** 2 * self.Tc ** 2 * a / self.Pc
        return a
    
    def PR_b(self):
        "计算PR方程的常数b"
        b = 0.077796 * R * self.Tc / self.Pc
        return b
    
    def PR_V(self, P, V, T, a, b):
        "求PR方程的体积根"
        V = b + R * T / ( P + a / ( V ** 2 + 2 * b * V - b ** 2 ) ) 
        return V
    
    def PR_PHI(self, P, V, T, a, b, Z):
        "求PR方程的phi"
        phi = Z - 1 - np.log( P * ( V - b ) / R / T) - \
        a / (2 ** 1.5) / b / R / T * \
        np.log( ( V + (2 ** 0.5 + 1) * b) / ( V - (2 ** 0.5 - 1) * b ))
        return phi
    
    def PR_PHI_MIX(self, P, V, T, a, b, xya, bi, Z):
        "求PR方程的 phi 混合物"
        phi = bi / b * (Z - 1) - np.log ( P * (V - b) / R / T) - a / ( 2 ** 1.5 ) / b / R / T * \
        ( 2 / a * xya - bi / b ) * np.log ( ( V + ( 2 ** 0.5 + 1) * b) / ( V - (2 ** 0.5 - 1) * b ) )
        phi = np.exp(phi).flatten()
        return phi

    def RK_a(self):
        "计算RK方程的常数a"
        self.RKa = 0.42748 * R ** 2 * self.Tc ** 2.5 / self.Pc
        return self.RKa
    
    def RK_b(self):
        "计算RK方程的常数b"
        self.RKb = 0.08664 * R * self.Tc / self.Pc
        return self.RKb

    def RK_V(self, P, V, T, a, b):
        "求RK方程的体积根"
        V = b + R * T / ( P + a / ( T ** 0.5 * V * ( V + b ) ) ) 
        return V